Show/hide code
# Parâmetros da distribuição Beta
a, b = 10, 4
median = round(beta.median(a, b), 3)
print(f"A mediana da distribuição Beta({a}, {b}): {median}")A mediana da distribuição Beta(10, 4): 0.725
O seu professor chega na sala de aula e mostra uma moeda. Você suspeita que a moeda possa ser falsa e ter duas caras. Considere a priori probabilidades iguais para os eventos da moeda ser falsa ou ser honesta (i.e. uma moeda bem equilibrada).
Para calcular a probabilidade de obter “cara” em um lançamento dessa moeda, considerando que a moeda pode ser honesta ou falsa. Suponha que:
Assumindo pelo enunciado:
A probabilidade de obter cara será dada por:
\(P(C_{1}) = P(C_{1}|H)P(H) + P(C_{1}|F)P(F)\)
Substituindo os valores fornecidos:
\(P(C_{1}) = (0.5 \times 0.5) + (1 \times 0.5)\)
\(P(C_{1}) = 0.25 + 0.5\)
\(P(C_{1}) = 0.75\)
Portanto, a probabilidade de obter cara em um lançamento dessa moeda é \(0.75\) ou 75%.
Queremos encontrar \(P(F|C_{1})\), a probabilidade de a moeda ser falsa dado que o resultado foi cara. Segue via Teorema de Bayes:
\(P(F|C_{1}) = \frac{P(C_{1}|F)P(F)}{P(C_{1})}\)
Logo:
\(P(F|C_{1}) = \frac{P(C_{1}|F)P(F)}{P(C_{1})} = \frac{1 \times 0.5}{0.75} = \frac{0.5}{0.75} = \frac{2}{3}\)
Portanto, a probabilidade de a moeda ser falsa, dado que o resultado do lançamento foi cara, é \(\frac{2}{3}\) ou aproximadamente 66.67%.
Assumimos:
Queremos calcular \(P(F|C^n)\), a probabilidade de a moeda ser falsa dado que todos os \(n\) lançamentos resultaram em cara.
Segue via Teorema de Bayes:
\(P(F|C^n) = \frac{P(C^n|F)P(F)}{P(C^n)}\)
Onde:
Para encontrar \(P(C^n)\), recorremos a probabilidade total:
\(P(C^n) = P(C^n|H)P(H) + P(C^n|F)P(F)\)
Substituindo os valores fornecidos:
\(P(C^n) = (0.5^n \times 0.5) + (1 \times 0.5)\)
\(P(C^n) = 0.5^{n+1} + 0.5\)
Agora, substituímos na fórmula de Bayes:
\(P(F|C^n) = \frac{1 \times 0.5}{0.5^{n+1} + 0.5}\)
$P(F|C^n) = # Inferência Bayesiana {.unnumbered}
Observa-se que valores grandes de \(n\), \(0.5^n\) se aproxima de 0 muito rapidamente, porque \(0.5^n\) decresce exponencialmente.
Portanto, para \(n\) grande, \(P(F|C^n)\) se aproxima de:
\(P(F|C^n) \approx \frac{1}{0 + 1} = 1\)
Isso significa que, conforme \(n\) aumenta, a probabilidade de a moeda ser falsa se aproxima de 1 (ou 100%).
Portanto, se o professor lançar a moeda \(n\) vezes e obter \(n\) caras, a probabilidade de a moeda ser falsa tende a 1 conforme \(n\) se torna grande.
Após um lançamento resultando em “cara”, a probabilidade de a moeda ser falsa, \(P(F|C)\), foi calculada como \(\frac{2}{3}\). Consequentemente, a probabilidade de a moeda ser honesta, \(P(H|C)\), será:
\(P(H|C) = 1 - P(F|C) = 1 - \frac{2}{3} = \frac{1}{3}\)
Para a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento foi “cara” basta obter a soma das probabilidades condicionais de obter “cara” no próximo lançamento, dadas as duas possibilidades sobre a moeda (falsa ou honesta):
\(P(C_2|C_1) = P(C_2|F)P(F|C_1) + P(C_2|H)P(H|C_1)\)
Onde:
Substituindo os valores, temos:
\(P(C_2|C_1) = (1 \times \frac{2}{3}) + (0.5 \times \frac{1}{3})\)
\(P(C_2|C_1) = \frac{2}{3} + \frac{1}{6}\)
\(P(C_2|C_1) = \frac{4}{6} + \frac{1}{6}\)
\(P(C_2|C_1) = \frac{5}{6}\)
Portanto, a probabilidade de obter “cara” no próximo lançamento, dado que o primeiro lançamento resultou em “cara”, é \(\frac{5}{6}\) ou aproximadamente 83.33%.
No contexto descrito, a afirmação “os dois lançamentos da moeda são independentes” é falsa porque o resultado de cada lançamento influencia nossa crença sobre a natureza da moeda (se ela é honesta ou falsa). Essa influência altera a probabilidade dos resultados subsequentes.
Vale ressaltar que dois eventos \(A\) e \(B\) são independentes se e somente se:
\(P(A \cap B) = P(A) \cdot P(B)\)
Ou, de forma equivalente:
\(P(B|A) = P(B)\)
Se os lançamentos fossem verdadeiramente independentes, o resultado de um lançamento não alteraria a probabilidade do resultado do próximo lançamento. No entanto, neste contexto, o resultado do primeiro lançamento afeta nossa crença sobre a natureza da moeda, que por sua vez afeta a probabilidade do segundo lançamento. Ou seja , os lançamentos são dependentes porque \(P(C_2 | C_1) \neq P(C_2)\). Aqui, \(C_1\) e \(C_2\) representam o primeiro e o segundo lançamento resultando em cara, respectivamente.
A afirmação correta seria: “Os dois lançamentos da moeda são condicionalmente dependentes”.
Em resumo, o resultado de um lançamento influencia a probabilidade do próximo lançamento devido à incerteza sobre a natureza da moeda. Isso demonstra que os lançamentos não são independentes; são, na verdade, condicionalmente dependentes.
Seja \(y_{1} , y_{2} ,\dots, y_{n}\) uma amostra da distribuição de Bernoulli com probabilidade de sucesso θ e considere uma distribuição a priori uniforme para θ.
Para uma amostra de Bernoulli, a verossimilhança é dada por:
\(P(y_1, y_2, \ldots, y_n | \theta) = \theta^{\sum_{i=1}^n y_i} (1 - \theta)^{n - \sum_{i=1}^n y_i}\)
Seja \(k = \sum_{i=1}^n y_i\), que é o número de sucessos na amostra. Assim, a função de verossimilhança pode ser escrita como:
\(P(y_1, y_2, \ldots, y_n | \theta) = \theta^k (1 - \theta)^{n - k}\)
A distribuição a priori uniforme para \(\theta\) no intervalo \([0, 1]\) é dada por:
\(P(\theta) = 1 \quad \text{para} \quad \theta \in [0, 1]\)
Distribuição a posteriori:
\(P(\theta | y_1, y_2, \ldots, y_n) \propto P(y_1, y_2, \ldots, y_n | \theta) \cdot P(\theta)\)
Substituindo as funções:
\(P(\theta | y_1, y_2, \ldots, y_n) \propto \theta^k (1 - \theta)^{n - k} \cdot 1\)
\(P(\theta | y_1, y_2, \ldots, y_n) \propto \theta^k (1 - \theta)^{n - k}\)
Essa é a forma da distribuição Beta. Portanto, a distribuição a posteriori de \(\theta\) é uma distribuição Beta com parâmetros \(\alpha = k + 1\) e \(\beta = n - k + 1\):
\(\theta | y_1, y_2, \ldots, y_n \sim \text{Beta}(k + 1, n - k + 1)\)
Para uma distribuição \(\text{Beta}(\alpha, \beta)\), a média e a variância são dadas por:
\(\mathbb{E}[\theta] = \frac{\alpha}{\alpha + \beta}\)
\(\text{Var}(\theta) = \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}\)
Substituindo os parâmetros \(\alpha = k + 1\) e \(\beta = n - k + 1\):
Média:
\(\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{(k + 1) + (n - k + 1)} = \frac{k + 1}{n + 2}\)
Variância:
\(\text{Var}(\theta | y_1, y_2, \ldots, y_n) = \frac{(k + 1)(n - k + 1)}{((k + 1) + (n - k + 1))^2 ((k + 1) + (n - k + 1) + 1)}\)
\(\text{Var}(\theta | y_1, y_2, \ldots, y_n) = \frac{(k + 1)(n - k + 1)}{(n + 2)^2 (n + 3)}\)
A distribuição a posteriori de \(\theta\) foi encontrada como \(\text{Beta}(k + 1, n - k + 1)\), tendo \(\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{n + 2}\).
A priori, temos uma distribuição uniforme para \(\theta\) no intervalo \([0, 1]\), cuja esperança é:
\(\mathbb{E}[\theta] = \frac{1}{2}\)
Sabendo que a estimativa de máxima verossimilhança (MLE) de \(\theta\) para uma amostra Bernoulli é a proporção de sucessos na amostra:
\(\hat{\theta} = \frac{k}{n}\)
Basta substituir os valores encontrados:
\(\mathbb{E}[\theta | y_1, y_2, \ldots, y_n] = \frac{k + 1}{n + 2}\)
\(\mathbb{E}[\theta] = \frac{1}{2}\)
\(\hat{\theta} = \frac{k}{n}\)
Vamos expressar \(\frac{k + 1}{n + 2}\) na forma desejada:
\(\frac{k + 1}{n + 2} = \left(1 - \frac{2}{n + 2}\right) \times \frac{1}{2} + \frac{2}{n + 2} \times \frac{k}{n}\)
Para verificar, expandimos a expressão:
\(\left(1 - \frac{2}{n + 2}\right) \times \frac{1}{2} + \frac{2}{n + 2} \times \frac{k}{n}\)
\(= \frac{1}{2} - \frac{1}{n + 2} + \frac{2k}{n(n + 2)}\)
\(\frac{k + 1}{n + 2} = \frac{n}{2(n + 2)} + \frac{2k}{n(n + 2)}\)
Efetuando mais uma manipulação algébrica:
\(\frac{k + 1}{n + 2} = \frac{1}{2} \left( \frac{n}{n + 2} \right) + \left( \frac{2}{n + 2} \right) \times \frac{k}{n}\)
Logo, \(w = \frac{2}{n + 2}\) e \((1 - w) = \frac{n}{n + 2}\). Isso significa que a esperança a posteriori de \(\theta\) é uma média ponderada entre a esperança a priori e a estimativa de máxima verossimilhança, com os pesos dependendo do tamanho da amostra \(n\).
Ou seja:
Quando a amostra é pequena (\(n\) pequeno), a esperança a priori tem um peso maior na esperança a posteriori.
Quando a amostra é grande (\(n\) grande), a estimativa de máxima verossimilhança tem um peso maior, e a esperança a posteriori se aproxima da MLE.
Isso reflete a intuição de que com mais dados, nossa crença sobre \(\theta\) é mais fortemente influenciada pela evidência observada do que pelas crenças a priori.
Distribuição preditiva de \(y_{n+1}\):
Para obter a distribuição preditiva de \(y_{n+1}\), precisamos calcular a probabilidade de \(y_{n+1}\) condicionado aos dados observados \(y_1, y_2, \dots, y_n\). A distribuição de \(y_{n+1}\) dado \(\theta\) é Bernoulli com parâmetro \(\theta\):
\(p(y_{n+1} | \theta) = \theta^{y_{n+1}} (1 - \theta)^{1 - y_{n+1}}\)
A distribuição preditiva é obtida integrando \(\theta\) fora desta expressão, ponderando pela distribuição a posteriori de \(\theta\), encontrada previamente:
\(p(y_{n+1} | y_1, y_2, \dots, y_n) = \int_0^1 p(y_{n+1} | \theta) p(\theta | y_1, y_2, \dots, y_n) \, d\theta\)
Substituindo \(p(y_{n+1} | \theta)\) e \(p(\theta | y_1, y_2, \dots, y_n)\):
\(p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \int_0^1 \theta \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)
\(p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \int_0^1 (1 - \theta) \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)
Para \(y_{n+1} = 1\):
\(p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \int_0^1 \theta \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)
\(= \frac{1}{B(k + 1, n - k + 1)} \int_0^1 \theta^{k + 1 - 1} (1 - \theta)^{n - k} \, d\theta\)
\(= \frac{B(k + 2, n - k + 1)}{B(k + 1, n - k + 1)}\)
\(= \frac{\frac{\Gamma(k + 2) \Gamma(n - k + 1)}{\Gamma(n + 3)}}{\frac{\Gamma(k + 1) \Gamma(n - k + 1)}{\Gamma(n + 2)}}\)
\(= \frac{(k + 1) \Gamma(k + 1) \Gamma(n - k + 1) / \Gamma(n + 3)}{\Gamma(k + 1) \Gamma(n - k + 1) / \Gamma(n + 2)}\)
\(= \frac{(k + 1)}{n + 2}\)
Para \(y_{n+1} = 0\):
\(p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \int_0^1 (1 - \theta) \cdot \frac{\theta^k (1 - \theta)^{n - k}}{B(k + 1, n - k + 1)} \, d\theta\)
\(= \frac{1}{B(k + 1, n - k + 1)} \int_0^1 \theta^k (1 - \theta)^{n - k + 1 - 1} \, d\theta\)
\(= \frac{B(k + 1, n - k + 2)}{B(k + 1, n - k + 1)}\)
\(= \frac{\frac{\Gamma(k + 1) \Gamma(n - k + 2)}{\Gamma(n + 3)}}{\frac{\Gamma(k + 1) \Gamma(n - k + 1)}{\Gamma(n + 2)}}\)
\(= \frac{\Gamma(n - k + 2) / \Gamma(n + 3)}{\Gamma(n - k + 1) / \Gamma(n + 2)}\)
\(= \frac{(n - k + 1)}{n + 2}\)
Assim, a distribuição preditiva \(p(y_{n+1} | y_1, y_2, \dots, y_n)\) é:
\(p(y_{n+1} = 1 | y_1, y_2, \dots, y_n) = \frac{k + 1}{n + 2}\)
\(p(y_{n+1} = 0 | y_1, y_2, \dots, y_n) = \frac{n - k + 1}{n + 2}\)
A distribuição preditiva de uma observação futura \(y_{n+1}\) é a média ponderada das proporções de sucessos e fracassos observadas, ajustadas pelo tamanho da amostra, refletindo o princípio de suavização de Laplace, onde mesmo com uma grande quantidade de dados, há um pequeno ajuste baseado na crença a priori.
No exercício 2, calcule:
Para calcular a estimativa Bayesiana sob a perda quadrática, precisamos encontrar o estimador que minimiza o valor esperado da perda quadrática, geralmente expressa como \(L(\theta, d) = (\theta - d)^2\), onde \(\theta\) é o parâmetro verdadeiro e \(d\) é a decisão ou estimativa tomada. Logo:
\(\mathbb{E}[(\theta - d)^2 | y_1, y_2, \dots, y_n]\)
onde \(\theta\) segue a distribuição a posteriori dado os dados observados \(y_1, y_2, \dots, y_n\).
Dado que \(\theta\) tem distribuição a posteriori \(\text{Beta}(k + 1, n - k + 1)\), a estimativa Bayesiana que minimiza a perda quadrática é a média da distribuição a posteriori, afinal, a média de uma distribuição é o estimador de mínimo erro quadrático médio (MSE).
A média da distribuição \(\text{Beta}(\alpha, \beta)\) é \(\frac{\alpha}{\alpha + \beta}\). Para a nossa distribuição a posteriori \(\text{Beta}(k + 1, n - k + 1)\):
\(\mathbb{E}[\theta | y_1, y_2, \dots, y_n] = \frac{k + 1}{(k + 1) + (n - k + 1)} = \frac{k + 1}{n + 2}\)
Portanto, a estimativa Bayesiana de \(\theta\) para a perda quadrática, dada a amostra de dados \([y_1, y_2, \dots, y_n]\), é:
\(d^* = \frac{k + 1}{n + 2}\)
A estimativa Bayesiana \(d^*\) corresponde à média a posteriori do parâmetro \(\theta\), levando em conta tanto a informação a priori (que \(\theta\) é uniformemente distribuído entre 0 e 1) quanto a informação proveniente da amostra (\(k\) sucessos em \(n\) tentativas). O ajuste de \(+1\) no numerador e \(+2\) no denominador é uma forma de regularização que evita extremos baseados em amostras pequenas, fornecendo uma suavização estatística que protege contra conclusões extremas com dados limitados, refletindo um compromisso entre o conhecimento a priori sobre a distribuição de \(\theta\) e a evidência obtida da amostra, alinhando-se com os princípios fundamentais da inferência Bayesiana.
A função de perda zero-um é dada por \(L(\theta, d) = 1\) se \(\theta \neq d\) e \(L(\theta, d) = 0\) se \(\theta = d\). Essencialmente, essa perda penaliza decisões incorretas sem considerar quão “erradas” elas estão, ao contrário da perda quadrática que penaliza de acordo com o quadrado da diferença entre a verdadeira \(\theta\) e a estimativa \(d\).
Cabe ressaltar que na estimativa para perda zero-um, procuramos minimizar o erro de classificação. Em geral, o estimador que minimiza a esperança da perda zero-um sob uma distribuição de probabilidades é o valor de \(\theta\) que maximiza a distribuição a posteriori, ou seja, a moda da distribuição a posteriori.
No nosso caso especial, a distribuição a posteriori de \(\theta\) após observar os dados é a \(\text{Beta}(s + 1, n - s + 1)\), onde \(s = \sum_{i=1}^n y_i\). Para \(n = 12\) e \(s = 9\), temos:
\(\theta | y_1, y_2, \dots, y_{12} \sim \text{Beta}(9 + 1, 12 - 9 + 1) = \text{Beta}(10, 4)\)
A moda de uma distribuição \(\text{Beta}(\alpha, \beta)\), para \(\alpha > 1\) e \(\beta > 1\), é dada por:
\(\text{Mode}[\theta] = \frac{\alpha - 1}{\alpha + \beta - 2}\)
Para a nossa distribuição \(\text{Beta}(10, 4)\):
\(\text{Mode}[\theta] = \frac{10 - 1}{10 + 4 - 2} = \frac{9}{12} = 0.75\)
O limite de \(\epsilon \rightarrow 0\) na pergunta sugere considerar um caso extremo da função de perda ou dos parâmetros envolvidos. Neste contexto, normalmente analisamos o comportamento da estimativa à medida que o tamanho da amostra cresce ou como a estimativa se comporta sob condições extremas da amostra ou dos parâmetros, todavia, como já calculamos a moda, o que seria o resultado para a perda zero-um, não há dependência direta de \(\epsilon\) no cálculo, a menos que \(\epsilon\) esteja relacionado a algum tipo de regularização ou parâmetro de suavização não especificado explicitamente na questão.
Para encontrar a estimativa Bayesiana sob a perda absoluta, precisamos minimizar o valor esperado da perda absoluta. A função de perda absoluta é dada por \(L(\theta, d) = |\theta - d|\), onde \(\theta\) é o parâmetro verdadeiro e \(d\) é a decisão ou estimativa tomada.
O estimador que minimiza a perda absoluta é a mediana da distribuição a posteriori, afinal, a mediana minimiza a soma das distâncias absolutas em relação a ela mesma.
A distribuição a posteriori de \(\theta\) após observar os dados $ y_1, y_2, , y_n $ é \(\text{Beta}(k + 1, n - k + 1)\), onde \(k = \sum_{i=1}^n y_i\). No caso especial onde $ n = 12 $ e $ k = s = 9 $, temos:
\(\theta | y_1, y_2, \ldots, y_{12} \sim \text{Beta}(10, 4)\)
Para uma distribuição \(\text{Beta}(\alpha, \beta)\), a mediana não tem uma fórmula fechada simples, mas pode ser aproximada ou calculada numericamente.
Uma aproximação comum para a mediana de uma distribuição \(\text{Beta}(\alpha, \beta)\) é:
\(\text{Mediana} \approx \frac{\alpha - \frac{1}{3}}{\alpha + \beta - \frac{2}{3}}\)
Substituindo \(\alpha = 10\) e \(\beta = 4\):
\(\text{Mediana} \approx \frac{10 - \frac{1}{3}}{10 + 4 - \frac{2}{3}} = \frac{9.6667}{13.3333} \approx 0.725\)
Portanto, a estimativa Bayesiana de \(\theta\) sob a perda absoluta, dada a amostra específica, é aproximadamente 0.725.
Ao optar por uma aproximação via python, obtemos:
# Parâmetros da distribuição Beta
a, b = 10, 4
median = round(beta.median(a, b), 3)
print(f"A mediana da distribuição Beta({a}, {b}): {median}")A mediana da distribuição Beta(10, 4): 0.725
De modo a encontrar o mesmo resultado para a estimativa Bayesiana de \(\theta\) sob a perda absoluta.
O intervalo de alta densidade posterior (HPD - Highest Posterior Density) é um intervalo que contém uma determinada percentagem da densidade de probabilidade e onde todas as probabilidades dentro do intervalo são maiores ou iguais às probabilidades fora do intervalo. Para calcular um intervalo HPD de 99% para a distribuição a posteriori de \(\theta\), onde \(\theta\) segue uma distribuição \(\text{Beta}(10, 4)\), utilizaremos o software python, pois a distribuição Beta não permite uma solução analítica simples para o intervalo HPD.
Para \(n = 12\) e \(s = 9\), a distribuição a posteriori de \(\theta\) é \(\text{Beta}(10, 4)\). A fórmula da função de densidade para a distribuição Beta é:
\(f(\theta; \alpha, \beta) = \frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)}\)
onde \(B(\alpha, \beta)\) é a função Beta, e para nosso caso específico, \(\alpha = 10\) e \(\beta = 4\).
Logo:
# Parâmetros da Beta
a, b = 10, 4
def hpd_interval(a, b, level=0.99):
lower_bound = (1 - level) / 2
upper_bound = 1 - lower_bound
lower_quantile = beta.ppf(lower_bound, a, b)
upper_quantile = beta.ppf(upper_bound, a, b)
return lower_quantile, upper_quantile
hpd_99 = hpd_interval(a, b, 0.99)
hpd_99_formatted = (round(hpd_99[0], 3), round(hpd_99[1], 3))
print(f"Intervalo HPD de 99% para Beta({a}, {b}): {hpd_99_formatted}")Intervalo HPD de 99% para Beta(10, 4): (0.379, 0.943)
O intervalo resultante explicita a região onde a densidade da distribuição Beta é mais alta e contém 99% da probabilidade total.
Suponha que \((x_{1} , x_{2} , x_{3})\) dado \(p_{1} , p_{2} , p_{3}\) segue uma distribuição Multinomial com parâmetros \(n\) e \((p_{1} , p_{2} , p_{3})\) , onde \(p_{i} ≥ 0 \space \text{e} \space p_{1} + p_{2} + p_{3} = 1\), e que, a priori, \((p_{1} , p_{2} , p_{3})\) segue uma distribuição de Dirichlet com parâmetros \((α_{1} , α_{2} , α_{3} )\):
A distribuição a posteriori de \((p_1, p_2, p_3)\) é uma distribuição de Dirichlet com parâmetros atualizados.
A distribuição Dirichlet a priori é:
\((p_1, p_2, p_3) \sim \text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3)\)
A função de verossimilhança para a distribuição multinomial é:
\(P(x_1, x_2, x_3 | p_1, p_2, p_3) = \frac{n!}{x_1! x_2! x_3!} p_1^{x_1} p_2^{x_2} p_3^{x_3}\)
Afinal, a probabilidade de obter exatamente \(x_1\) ocorrências do resultado 1, \(x_2\) ocorrências do resultado 2, e \(x_3\) ocorrências do resultado 3 é dada pelo produto das probabilidades individuais elevadas ao número de ocorrências, multiplicado pelo número de permutações possíveis.
A distribuição a posteriori será:
\(P(p_1, p_2, p_3 | x_1, x_2, x_3) \propto p_1^{\alpha_1 - 1} p_2^{\alpha_2 - 1} p_3^{\alpha_3 - 1} \cdot p_1^{x_1} p_2^{x_2} p_3^{x_3}\)
\(\propto p_1^{\alpha_1 + x_1 - 1} p_2^{\alpha_2 + x_2 - 1} p_3^{\alpha_3 + x_3 - 1}\)
Portanto, a distribuição a posteriori é uma distribuição de Dirichlet com parâmetros \((\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)\):
\((p_1, p_2, p_3) | (x_1, x_2, x_3) \sim \text{Dirichlet}(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)\)
As distribuições marginais de \(p_i\) para uma distribuição de Dirichlet são Beta. Para \(p_1\):
\(p_1 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_1 + x_1, \alpha_2 + x_2 + \alpha_3 + x_3)\)
Similarmente, para \(p_2\):
\(p_2 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_2 + x_2, \alpha_1 + x_1 + \alpha_3 + x_3)\)
E para \(p_3\):
\(p_3 | (x_1, x_2, x_3) \sim \text{Beta}(\alpha_3 + x_3, \alpha_1 + x_1 + \alpha_2 + x_2)\)
Estimativas Bayesianas de \(p_i\):
Para a distribuição Dirichlet \(\text{Dirichlet}(\alpha_1, \alpha_2, \alpha_3)\), a esperança de \(p_i\) é:
\(\mathbb{E}[p_i | x_1, x_2, x_3] = \frac{\alpha_i + x_i}{\sum_{j=1}^3 (\alpha_j + x_j)}\)
Portanto, as estimativas Bayesianas de \(p_1\), \(p_2\), e \(p_3\) são:
\(\hat{p_1} = \frac{\alpha_1 + x_1}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\) \(\hat{p_2} = \frac{\alpha_2 + x_2}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\) \(\hat{p_3} = \frac{\alpha_3 + x_3}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\)
Estimativa de \(p_j - p_i\):
Para calcular a estimativa Bayesiana de \(p_j - p_i\) sob a perda quadrática, usaremos a linearidade da esperança. A estimativa é simplesmente a diferença das esperanças de \(p_j\) e \(p_i\):
\(\mathbb{E}[p_j - p_i | x_1, x_2, x_3] = \mathbb{E}[p_j | x_1, x_2, x_3] - \mathbb{E}[p_i | x_1, x_2, x_3]\)
Portanto:
\(\mathbb{E}[p_j - p_i | x_1, x_2, x_3] = \frac{\alpha_j + x_j}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3} - \frac{\alpha_i + x_i}{\alpha_1 + x_1 + \alpha_2 + x_2 + \alpha_3 + x_3}\)
Podemos assim, concatenar os resultados obtidos :
A distribuição a posteriori de \((p_1, p_2, p_3)\) é \(\text{Dirichlet}(\alpha_1 + x_1, \alpha_2 + x_2, \alpha_3 + x_3)\).
As distribuições marginais de \(p_i\) são \(\text{Beta}(\alpha_i + x_i, \sum_{j \neq i} (\alpha_j + x_j))\).
As estimativas Bayesianas de \(p_i\) são \(\hat{p_i} = \frac{\alpha_i + x_i}{\sum_{j=1}^3 (\alpha_j + x_j)}\).
As estimativas Bayesianas de \(p_j - p_i\) são \(\hat{p_j} - \hat{p_i} = \frac{\alpha_j + x_j}{\sum_{k=1}^3 (\alpha_k + x_k)} - \frac{\alpha_i + x_i}{\sum_{k=1}^3 (\alpha_k + x_k)}\).
É conhecido que 25% dos pacientes de um certo grupo que sofrem de enxaqueca melhoram após duas horas de serem tratados com um placebo. Para verificar se uma droga nova é melhor que o placebo, \(n = 20\) pacientes foram tratados com o placebo e verificou-se que após duas horas \(s = 8\) deles relataram ter melhorado. Seja \(\theta\) a probabilidade de um paciente tratado com a droga nova melhorar após duas horas.
A distribuição a priori uniforme \(\theta \sim \text{Uniforme}(0, 1)\):
\(P(\theta \leq 0.25) = 0.25\)
\(P(\theta > 0.25) = 1 - 0.25 = 0.75\)
As chances relativas a priori:
\(O(H_1) = \frac{P(\theta > 0.25)}{P(\theta \leq 0.25)} = \frac{0.75}{0.25} = 3\)
Sabemos que a distribuição a posteriori é \(\theta | s \sim \text{Beta}(9, 13)\).
Probabilidade a Posteriori de $ > 0.25 $:
\(P(\theta > 0.25 | s) = 1 - F_{\text{Beta}}(0.25 | 9, 13)\)
Consequentemente:
\(P(\theta \leq 0.25 | s) = F_{\text{Beta}}(0.25 | 9, 13)\)
Via python:
# Parâmetros da distribuição Beta
a, b = 9, 13
posterior_prob_H1 = 1 - beta.cdf(0.25, a, b)
posterior_prob_H0 = beta.cdf(0.25, a, b)
print(f"P(θ > 0.25 | s) = {posterior_prob_H1:.3f}")
print(f"P(θ ≤ 0.25 | s) = {posterior_prob_H0:.3f}")P(θ > 0.25 | s) = 0.944
P(θ ≤ 0.25 | s) = 0.056
Chances Relativas a Posteriori:
\(O(H_1 | s) = \frac{P(\theta > 0.25 | s)}{P(\theta \leq 0.25 | s)} = \frac{0.944}{0.056} = 16.810\)
Fator de Bayes:
\(B_{10}(s) = \frac{O(H_1 | s)}{O(H_1)} = \frac{16.810}{3} = 5.603\)
Definindo a função de perda de Neyman:
Considerando a função de perda de Neyman:
A decisão ótima é rejeitar \(H_0\) se:
\(P(\theta > 0.25 | s) > \frac{a_0}{a_0 + a_1} = \frac{5a_1}{5a_1 + a_1} = \frac{5}{6} = 0.833\)
Como \(P(\theta > 0.25 | s) = 0.944\) é maior do que \(0.833\), a decisão ótima é rejeitar \(H_0\).
Uma distribuição uniforme pode ser vista como “não informativa” em um certo sentido, mas pode não ser a mais adequada. Para hipóteses específicas como \(H_0\) e \(H_1\), uma distribuição Beta poderia ser mais informativa.
A priori uniforme favorece \(H_1\) no sentido que \(P(\theta > 0.25) = 3P(\theta \leq 0.25)\). Uma forma de corrigir isso é usar uma distribuição Beta com \(\beta = 3\alpha\), o que dá lugar a uma priori imprópria alternativa:
\(\theta \sim \text{Beta}(\alpha, 3\alpha)\)
Para obter uma variância maximizada, precisamos definir \(\alpha \to 0\). Então:
\(P(\theta) \propto \frac{1}{\theta(1 - \theta)}\)
Com a priori imprópria \(P(\theta) \propto \frac{1}{\theta(1 - \theta)}\), a distribuição a posteriori é proporcional a:
\(P(\theta | s) \propto P(s | \theta) P(\theta) \propto \theta^s (1 - \theta)^{n - s} \cdot \frac{1}{\theta(1 - \theta)} = \theta^{s-1} (1 - \theta)^{n-s-1}\)
Isto é uma distribuição \(\text{Beta}(s, n-s)\).
Para \(n = 20\) e \(s = 8\):
\(\theta | s \sim \text{Beta}(8, 12)\)
alpha_prior, beta_prior = 1, 3
n = 20
s = 8
a_post = alpha_prior + s
b_post = beta_prior + n - s
posterior_prob_H1_beta = 1 - beta.cdf(0.25, a_post, b_post)
print(f"P(θ > 0.25 | s) = {posterior_prob_H1_beta:.3f}")P(θ > 0.25 | s) = 0.904
A decisão de rejeitar \(H_0\) é mantida tanto com a priori uniforme quanto com a priori Beta(1, 3).
Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição de Poisson com média \(\theta\), e considere a priori que \(\theta\) tem uma distribuição Gama com parâmetros \(\alpha\) e \(\beta\) (isto é, com média \(\alpha\beta^{-1}\) e variância \(\alpha\beta^{-2}\)).
Dada uma amostra \(x_1, x_2, \ldots, x_n\) da distribuição de Poisson com média \(\theta\), a verossimilhança é:
\(P(x_1, x_2, \ldots, x_n | \theta) = \prod_{i=1}^n \frac{\theta^{x_i} e^{-\theta}}{x_i!}\)
\(P(x_1, x_2, \ldots, x_n | \theta) = \theta^{\sum_{i=1}^n x_i} e^{-n\theta} \prod_{i=1}^n \frac{1}{x_i!}\)
Considerando que \(\theta\) tem uma distribuição Gama a priori com parâmetros \(\alpha\) e \(\beta\):
\(P(\theta) \propto \theta^{\alpha - 1} e^{-\beta \theta}\)
A distribuição a posteriori é:
\(P(\theta | x_1, x_2, \ldots, x_n) \propto P(x_1, x_2, \ldots, x_n | \theta) P(\theta)\)
Substituindo as funções de verossimilhança e a priori:
\(P(\theta | x_1, x_2, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i} e^{-n\theta} \cdot \theta^{\alpha - 1} e^{-\beta \theta}\)
\(P(\theta | x_1, x_2, \ldots, x_n) \propto \theta^{(\alpha + \sum_{i=1}^n x_i) - 1} e^{-(\beta + n) \theta}\)
Portanto, a distribuição a posteriori é uma distribuição Gama com parâmetros:
\(\alpha^* = \alpha + \sum_{i=1}^n x_i\)
\(\beta^* = \beta + n\)
A média e variância da distribuição a posteriori são:
\(\mathbb{E}[\theta | x] = \frac{\alpha^*}{\beta^*} = \frac{\alpha + \sum_{i=1}^n x_i}{\beta + n}\)
\(\text{Var}[\theta | x] = \frac{\alpha^*}{(\beta^*)^2} = \frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2}\)
Sabemos que:
\(\mathbb{E}[\theta | x] = \frac{\alpha + \sum_{i=1}^n x_i}{\beta + n}\)
Podemos escrever:
\(\mathbb{E}[\theta | x] = \frac{\alpha}{\beta + n} + \frac{\sum_{i=1}^n x_i}{\beta + n}\)
Definimos \(\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\), assim:
\(\mathbb{E}[\theta | x] = \frac{\alpha}{\beta + n} + \frac{n\bar{x}}{\beta + n}\)
Se definirmos \(w = \frac{n}{\beta + n}\), temos:
\(\mathbb{E}[\theta | x] = (1 - w) \frac{\alpha}{\beta} + w\bar{x}\)
Portanto, podemos escrever:
\(\mathbb{E}[\theta | x] = w\bar{x} + (1-w)\alpha\beta^{-1}\)
Isso mostra que a estimativa a posteriori é uma combinação ponderada entre a média a priori \(\alpha\beta^{-1}\) e a média amostral \(\bar{x}\).
Quando \(\beta\) é grande com \(\alpha\beta^{-1}\) fixo:
\(w = \frac{n}{\beta + n} \approx \frac{n}{\beta}\)
Se \(\beta \to \infty\), então \(w \to 0\), o que significa que a esperança a posteriori \(\mathbb{E}[\theta | x]\) será dominada pela média a priori \(\alpha\beta^{-1}\), resultando em:
\(\mathbb{E}[\theta | x] \approx (1 - 0) \alpha\beta^{-1} = \alpha\beta^{-1}\)
A variância a posteriori é:
\(\text{Var}[\theta | x] = \frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2}\)
A variância a priori é:
\(\text{Var}[\theta] = \frac{\alpha}{\beta^2}\)
Queremos encontrar \(c\) tal que \(\text{Var}[\theta | x] > \text{Var}[\theta]\) quando \(\bar{x} > c\).
\(\frac{\alpha + \sum_{i=1}^n x_i}{(\beta + n)^2} > \frac{\alpha}{\beta^2}\)
Multiplicando ambos os lados por \((\beta + n)^2 \beta^2\):
\((\alpha + \sum_{i=1}^n x_i) \beta^2 > \alpha (\beta + n)^2\)
\(\alpha \beta^2 + \sum_{i=1}^n x_i \beta^2 > \alpha (\beta^2 + 2\beta n + n^2)\)
\(\sum_{i=1}^n x_i \beta^2 > \alpha (2\beta n + n^2)\)
Dividindo ambos os lados por \(n \beta^2\):
\(\bar{x} > \frac{\alpha (2\beta n + n^2)}{n \beta^2}\)
\(\bar{x} > \frac{\alpha (2\beta + n)}{\beta^2}\)
Portanto, o valor crítico \(c\) é:
\(c = \frac{\alpha (2\beta + n)}{\beta^2} = \frac{\alpha}{\beta} \times (2 + \frac{n}{\beta}) = (2 + \frac{n}{\beta}) \times \mathbb{E}[\theta]\)
A estimativa sob perda quadrática é imediata: \(\mathbb{E}[\theta | D] = \alpha^* / \beta^* = (s + 1) / (n + 1)\).
O limite da estimativa sob perda zero-um é a moda da distribuição a posteriori, que nesse caso é \((\alpha^* - 1) / \beta^* = s / (n + 1)\).
Agora, quando \(n = 10\) e \(\bar{x} = 1.55\), temos que \(\alpha^* = 16.50\) e \(\beta^* = 11\), de forma que a estimativa sob PA é a mediana \(\text{med}(\theta | D) \approx 1.470\).
Para o caso \(n = 10\) e \(\bar{x} = 1.55\), ache:
Temos que \(\alpha^* = 16.50\) e \(\beta^* = 11\), de forma que a estimativa sob PA é a mediana \(\text{med}(\theta | D) \approx 1.470\).
Vamos implementar os cálculos em Python para verificar os resultados obtidos anteriormente e, por fim, explicitar o intervalo HPD:
n = 10
x_bar = 1.55
s = n * x_bar
alpha_post = 1 + s
beta_post = 1 + n
theta_quadratic = alpha_post / beta_post
print(f"Estimativa bayesiana sob perda quadrática: {theta_quadratic:.3f}")
theta_zero_one = (alpha_post - 1) / beta_post
print(f"Limite da estimativa sob perda zero-um: {theta_zero_one:.3f}")
theta_absolute = gamma.median(alpha_post, scale=1/beta_post)
print(f"Estimativa bayesiana sob perda absoluta: {theta_absolute:.3f}")
hpd_interval = gamma.interval(0.95, alpha_post, scale=1/beta_post)
print(f"Intervalo HPD de 95%: ({hpd_interval[0]:.3f}, {hpd_interval[1]:.3f})")Estimativa bayesiana sob perda quadrática: 1.500
Limite da estimativa sob perda zero-um: 1.409
Estimativa bayesiana sob perda absoluta: 1.470
Intervalo HPD de 95%: (0.866, 2.306)
Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição Normal com média \(\mu\) e variância \(\phi^{-1}\) conhecida, e considere a distribuição a priori \(\mu \sim \mathcal{N}(\mu_0, \tau^{-1})\).
Dada a amostra \(x_1, x_2, \ldots, x_n\) da distribuição Normal com média \(\mu\) e variância \(\phi^{-1}\), a verossimilhança é:
\(f(x_i | \mu) = \sqrt{\frac{\phi}{2\pi}} \exp\left( -\frac{\phi}{2}(x_i - \mu)^2 \right)\)
A função de verossimilhança conjunta para a amostra é:
\(L(\mu | x_1, x_2, \ldots, x_n) = \left(\frac{\phi}{2\pi}\right)^{n/2} \exp\left( -\frac{\phi}{2} \sum_{i=1}^n (x_i - \mu)^2 \right)\)
A priori \(\mu \sim \mathcal{N}(\mu_0, \tau^{-1})\) é dada por:
\(f(\mu) = \sqrt{\frac{\tau}{2\pi}} \exp\left( -\frac{\tau}{2} (\mu - \mu_0)^2 \right)\)
Onde teremos a distribuição a posteriori:
\(f(\mu | x_1, x_2, \ldots, x_n) \propto L(\mu | x_1, x_2, \ldots, x_n) f(\mu)\)
Substituindo as expressões:
\(f(\mu | x_1, x_2, \ldots, x_n) \propto \exp\left( -\frac{\phi}{2} \sum_{i=1}^n (x_i - \mu)^2 \right) \exp\left( -\frac{\tau}{2} (\mu - \mu_0)^2 \right)\)
Agrupando os termos quadráticos em \(\mu\):
\(\propto \exp\left( -\frac{1}{2} \left[ \phi \sum_{i=1}^n x_i^2 + \tau \mu_0^2 - 2\mu \left( \phi \sum_{i=1}^n x_i + \tau \mu_0 \right) + \mu^2 (\phi n + \tau) \right] \right)\)
Completando o quadrado:
\(f(\mu | x_1, x_2, \ldots, x_n) \propto \exp\left( -\frac{1}{2} \left[ (\phi n + \tau)\left(\mu - \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\right)^2 \right] \right)\)
Portanto, a distribuição a posteriori é uma Normal com os parâmetros atualizados:
\(\mu | x_1, x_2, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} \right)\)
A esperança a posteriori é:
\(\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\)
Podemos reescrever essa expressão como uma combinação ponderada entre a média amostral \(\bar{x}\) e a média a priori \(\mu_0\). Definimos \(w = \frac{\phi n}{\phi n + \tau}\):
\(\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = w \bar{x} + (1 - w) \mu_0\)
onde \(\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\).
Este resultado mostra que a estimativa a posteriori é uma média ponderada entre a média da amostra \(\bar{x}\) e a média a priori \(\mu_0\), com pesos \(w\) e \(1 - w\) respectivamente. O peso \(w\) depende do tamanho da amostra e dos parâmetros \(\phi\) e \(\tau\).
Se \(\bar{x}_m\) é a média de \(m\) observações futuras \(x_{n+1}, \ldots, x_{n+m}\), condicionalmente independentes de \(x_1, \ldots, x_n\), a distribuição preditiva é:
\(p(\bar{x}_m | x_1, \ldots, x_n)\)
A média de \(m\) observações futuras é normalmente distribuída com a mesma média \(\mu\) e variância \(\frac{1}{\phi m}\). Portanto, a média preditiva é uma média ponderada da distribuição posterior de \(\mu\) e da distribuição das observações futuras.
A variância total é a soma das variâncias individuais:
\(\text{Var}(\bar{x}_m | x_1, \ldots, x_n) = \text{Var}(\mu | x_1, \ldots, x_n) + \text{Var}(x_{n+1}, \ldots, x_{n+m} | \mu)\)
\(= \frac{1}{\phi n + \tau} + \frac{1}{\phi m}\)
A distribuição preditiva de \(\bar{x}_m\) é:
\(\bar{x}_m | x_1, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi n \bar{x} + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} + \frac{1}{\phi m} \right)\)
Quando \(\tau \to 0\), a priori se torna não informativa, ou seja, \(p(\mu) \propto 1\). Nesse caso, o peso \(w\) na expressão da esperança a posteriori se torna 1, pois:
\(w = \frac{\phi n}{\phi n + \tau} \approx 1\)
Portanto, a esperança a posteriori se torna a média amostral:
\(\mathbb{E}[\mu | x_1, x_2, \ldots, x_n] \approx \bar{x}\)
E a variância a posteriori se torna:
\(\text{Var}(\mu | x_1, x_2, \ldots, x_n) \approx \frac{1}{\phi n}\)
Isso mostra que, com uma priori não informativa, a distribuição a posteriori é dominada pelos dados da amostra, e a média a posteriori é simplesmente a média amostral.
Seja \(x_1, x_2, \ldots, x_n\) uma amostra da distribuição Normal com média \(\mu\) e variância \(\phi^{-1}\) conhecida, e considere a distribuição a priori \(\mu \sim \mathcal{N}(\mu_0, \tau^{-1})\).
(a) Ache o estimador bayesiano de \(\mu\) no caso de (i) perda quadrática \((L(\hat{\mu}, \mu) = (\hat{\mu} - \mu)^2)\), (ii) perda absoluta \((L(\hat{\mu}, \mu) = |\hat{\mu} - \mu|)\) e (iii) perda zero-um \((L(\hat{\mu}, \mu) = I(|\hat{\mu} - \mu| \geq \epsilon))\).
Para uma função de perda quadrática \(L(\hat{\mu}, \mu) = (\hat{\mu} - \mu)^2\), o estimador bayesiano é a média da distribuição a posteriori.
Como sabemos da questão anterior que a distribuição a posteriori é uma Normal com parâmetros atualizados:
\(\mu | x_1, x_2, \ldots, x_n \sim \mathcal{N}\left( \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}, \frac{1}{\phi n + \tau} \right)\)
A média dessa distribuição é:
\(\hat{\mu}_{\text{PQ}} = \mathbb{E}[\mu | x_1, x_2, \ldots, x_n] = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\)
Para a distribuição normal, a média, a mediana e a moda são todos iguais devido à sua simetria. Assim, os estimadores sob perda quadrática (i), perda absoluta (ii) e perda zero-um (iii) são todos iguais a:
\(\hat{\mu} = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\)
Essa equivalência é uma consequência direta da simetria da distribuição normal, onde a distribuição é perfeitamente simétrica em torno da média, que coincide com a mediana e a moda.
(b) Ache o intervalo HPD para \(\mu\) com nível \(100(1 - \alpha)\%\).
O intervalo HPD (Highest Posterior Density) para uma distribuição Normal \(\mathcal{N}(\mu_{\text{post}}, \sigma^2_{\text{post}})\) é encontrado nos quantis da distribuição que correspondem ao nível de confiança desejado.
Os limites do intervalo HPD de \(100(1 - \alpha)\%\) são:
\(\left[ \mu_{\text{post}} - z_{\alpha/2} \sigma_{\text{post}}, \mu_{\text{post}} + z_{\alpha/2} \sigma_{\text{post}} \right]\)
onde \(\mu_{\text{post}} = \frac{\phi \sum_{i=1}^n x_i + \tau \mu_0}{\phi n + \tau}\) e \(\sigma^2_{\text{post}} = \frac{1}{\phi n + \tau}\).
Vamos calcular esses valores:
n = 10
x_bar = 1.55
phi = 1
tau = 1
mu_0 = 0
sum_x = n * x_bar
mu_post = (phi * sum_x + tau * mu_0) / (phi * n + tau)
sigma_post = np.sqrt(1 / (phi * n + tau))
alpha = 0.05
z = norm.ppf(1 - alpha / 2)
hpd_interval = [mu_post - z * sigma_post, mu_post + z * sigma_post]
print(f"Intervalo HPD de 95%: [{hpd_interval[0]:.3f}, {hpd_interval[1]:.3f}]")Intervalo HPD de 95%: [0.818, 2.000]
Considere o modelo \(X \sim \text{Uniforme} (\theta, \theta + 1)\) \((-\infty < \theta < \infty)\).
(a) Mostre que \(\theta\) é um parâmetro de locação.
Um parâmetro \(\theta\) é um parâmetro de locação se, ao adicionar uma constante \(c\) a \(\theta\), a distribuição resultante for uma translação da distribuição original.
Para o modelo \(X \sim \text{Uniforme} (\theta, \theta + 1)\):
\[ f_X(x; \theta) = \begin{cases} 1 & \text{se } \theta \leq x \leq \theta + 1 \\ 0 & \text{caso contrário} \end{cases} \]
Se adicionarmos uma constante \(c\) a \(\theta\), a nova distribuição será \(X \sim \text{Uniforme} (\theta + c, \theta + c + 1)\):
\[ f_X(x; \theta + c) = \begin{cases} 1 & \text{se } \theta + c \leq x \leq \theta + c + 1 \\ 0 & \text{caso contrário} \end{cases} \]
Essa é uma translação da distribuição original, mostrando que \(\theta\) é um parâmetro de locação.
(b) Dada a priori não informativa usual para o modelo de locação, \(p(\theta) \propto 1\), e uma amostra aleatória \(X_1, \ldots, X_n\) do modelo acima, discuta a propriedade da distribuição a posteriori.
A priori não informativa usual para um modelo de locação é \(p(\theta) \propto 1\).
Dada uma amostra aleatória \(X_1, \ldots, X_n\) do modelo \(\text{Uniforme} (\theta, \theta + 1)\):
A verossimilhança é:
\(L(\theta | X_1, \ldots, X_n) \propto \prod_{i=1}^n f_X(x_i; \theta)\)
\[ L(\theta | X_1, \ldots, X_n) = \begin{cases} 1 & \text{se } \theta \leq \min(X_i) \text{ e } \max(X_i) \leq \theta + 1 \\ 0 & \text{caso contrário} \end{cases} \]
Dado que a priori é \(p(\theta) \propto 1\), a distribuição a posteriori é proporcional à verossimilhança:
\(p(\theta | X_1, \ldots, X_n) \propto L(\theta | X_1, \ldots, X_n)\)
\[ p(\theta | X_1, \ldots, X_n) = \begin{cases} \frac{1}{1 - (\max(X_i) - \min(X_i))} & \text{se } \max(X_i) - 1 \leq \theta \leq \min(X_i) \\ 0 & \text{caso contrário} \end{cases} \]
Portanto, a distribuição a posteriori de \(\theta\) é uma distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\).
(c) Na situação da parte (b), ache:
Dado que a distribuição a posteriori é uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\), a média é:
\(\hat{\theta}_{\text{PQ}} = \frac{(\max(X_i) - 1) + \min(X_i)}{2}\)
A mediana de uma distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\) é:
\(\hat{\theta}_{\text{PA}} = \frac{(\max(X_i) - 1) + \min(X_i)}{2}\)
Para a distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\), o intervalo HPD \(100(1 - \alpha)\%\) é simplesmente o intervalo central que cobre \(100(1 - \alpha)\%\) da distribuição:
\(\text{HPD} = \left[ X_n - 1 + \frac{\alpha}{2} (X_1 - (X_n - 1)), X_1 - \frac{\alpha}{2} (X_1 - (X_n - 1)) \right]\)
Onde \(\min(X_i) = X_1\) e \(\max(X_i) = X_n\)
Para encontrar a probabilidade a posteriori do evento \(\theta \geq \theta_0\):
\(P(\theta \geq \theta_0 | X_1, \ldots, X_n)\)
Para a distribuição uniforme no intervalo \([ \max(X_i) - 1, \min(X_i)]\):
Se \(\theta_0 < \max(X_i) - 1\), então $ P(_0 | X_1, , X_n) = 1$.
Se \(\theta_0 > \min(X_i)\), então $ P(_0 | X_1, , X_n) = 0$.
Se \(\max(X_i) - 1 \leq \theta_0 \leq \min(X_i)\):
\(P(\theta \geq \theta_0 | X_1, \ldots, X_n) = \frac{\min(X_i) - \theta_0}{\min(X_i) - (\max(X_i) - 1)}\)
Considere o modelo \(X \sim \text{Uniforme} (0, \theta)\) \((0 < \theta < \infty)\).
Considere \(X \sim \text{Uniforme}(0, \theta)\). A função densidade de probabilidade (pdf) de \(X\) é:
\(f_X(x; \theta) = \theta^{-1} \mathbb{I}_{\{0 \leq x \leq \theta\}}\)
Vamos redefinir a distribuição em termos de \(\sigma\). Como \(\sigma = \theta^{-1}\), então \(\theta = \sigma^{-1}\). A pdf em termos de \(\sigma\) é:
\(f_X(x; \sigma) = \sigma \mathbb{I}_{\{0 \leq x \leq \sigma^{-1}\}}\)
Vamos multiplicar \(\sigma\) por uma constante \(c\). A nova densidade será:
\(f_Y(y; c\sigma) = c\sigma \mathbb{I}_{\{0 \leq y \leq (c\sigma)^{-1}\}}\)
Essa é uma escala da densidade original, mostrando que \(\sigma = \theta^{-1}\) é um parâmetro de escala.
Segue que:
\(L(\theta | X_1, \ldots, X_n) = \prod_{i=1}^n f_X(x_i; \theta)\)
Como \(X_i \sim \text{Uniforme}(0, \theta)\), a verossimilhança é:
\(L(\theta | X_1, \ldots, X_n) = \theta^{-n} \mathbb{I}_{\{0 \leq \max(X_i) \leq \theta\}}\)
A priori não informativa é \(p(\sigma) \propto \sigma^{-1}\). Como \(\sigma = \theta^{-1}\), temos \(p(\theta) \propto \theta^{-1}\). Logo, a distribuição posteriori será:
\(p(\theta | X_1, \ldots, X_n) \propto L(\theta | X_1, \ldots, X_n) p(\theta)\)
\(p(\theta | X_1, \ldots, X_n) \propto \theta^{-n} \mathbb{I}_{\{0 \leq \max(X_i) \leq \theta\}} \theta^{-1}\)
\(p(\theta | X_1, \ldots, X_n) \propto \theta^{-(n+1)} \mathbb{I}_{\{\theta \geq \max(X_i)\}}\)
Por fim:
\(p(\theta | X_1, \ldots, X_n) = \frac{(n \max(X_i))^n}{\theta^{n+1}} \mathbb{I}_{\{\theta \geq \max(X_i)\}}\)
Esta é uma distribuição de Pareto com forma \(n\) e escala \(\max(X_i)\).
Na situação da parte (b), ache:
os estimadores bayesianos sob Perda Quadrática e sob Perda Absoluta.
Perda Quadrática: Para a função de perda quadrática \(L(\hat{\theta}, \theta) = (\hat{\theta} - \theta)^2\), o estimador bayesiano é a média da distribuição a posteriori. A média de uma distribuição Pareto \(\text{Pareto}(\max(X_i), n)\) é dada por: \(\hat{\theta}_{\text{PQ}} = \frac{n \max(X_i)}{n - 1}\) , para \(n > 1\).
Perda Absoluta: Para a função de perda absoluta \(L(\hat{\theta}, \theta) = |\hat{\theta} - \theta|\), o estimador bayesiano é a mediana da distribuição a posteriori. A mediana de uma distribuição Pareto \(\text{Pareto}(\max(X_i), n)\) é dada por:
\(\hat{\theta}_{\text{PA}} = \max(X_i) \cdot 2^{1/(n)}\)
Para uma distribuição Pareto \(\text{Pareto}(\max(X_i), n+1)\), o intervalo HPD \(100(1 - \alpha)\%\) pode ser calculado encontrando os quantis apropriados. Seja \(F_{\text{Pareto}}^{-1}(p)\) a função quantil da distribuição Pareto:
\(\text{HPD} = [ \max(X_i), \max(X_i) (1 - \alpha)^{-\frac{1}{n+1}} ]\)
A probabilidade a posteriori do evento \(\theta > \theta_0\) é dada pela função de sobrevivência da distribuição Pareto:
\(P(\theta > \theta_0 | X_1, \ldots, X_n) = \left( \frac{\max(X_i)}{\theta_0} \right)^{n} \mathbb{I}_{\{\theta_0 \geq \max(X_i)\}}\)
Considere o modelo \(x_1, \ldots, x_n | \theta \sim \text{Ber}(\theta)\).
A priori de Jeffreys é uma escolha não informativa que é proporcional à raiz quadrada do determinante da matriz de informação de Fisher. Para o modelo Bernoulli \(\text{Ber}(\theta)\), vamos calcular a matriz de informação de Fisher e encontrar a priori de Jeffreys.
Para uma amostra \(x_1, x_2, \ldots, x_n\) de uma distribuição Bernoulli com parâmetro \(\theta\), a função de verossimilhança é:
\(L(\theta) = \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i}\)
A log-verossimilhança é:
\(\log L(\theta) = \left(\sum_{i=1}^n x_i\right) \log \theta + \left(n - \sum_{i=1}^n x_i\right) \log (1 - \theta)\)
A informação de Fisher é a expectativa do quadrado da derivada da log-verossimilhança em relação a \(\theta\):
\(I(\theta) = -\mathbb{E}\left[\frac{\partial^2}{\partial \theta^2} \log L(\theta)\right]\)
Calculando a segunda derivada da log-verossimilhança:
\(\frac{\partial}{\partial \theta} \log L(\theta) = \frac{\sum_{i=1}^n x_i}{\theta} - \frac{n - \sum_{i=1}^n x_i}{1 - \theta}\)
\(\frac{\partial^2}{\partial \theta^2} \log L(\theta) = -\frac{\sum_{i=1}^n x_i}{\theta^2} - \frac{n - \sum_{i=1}^n x_i}{(1 - \theta)^2}\)
Tomando a expectativa, considerando que \(\mathbb{E}[x_i] = \theta\):
\(I(\theta) = n \left(\frac{1}{\theta} + \frac{1}{1 - \theta}\right)\)
A priori de Jeffreys é proporcional à raiz quadrada da informação de Fisher:
\(\pi_J(\theta) \propto \sqrt{I(\theta)}\)
\(\pi_J(\theta) \propto \sqrt{n \left(\frac{1}{\theta} + \frac{1}{1 - \theta}\right)}\)
\(\pi_J(\theta) \propto \sqrt{\frac{n}{\theta(1 - \theta)}}\)
\(\pi_J(\theta) \propto \frac{1}{\sqrt{\theta(1 - \theta)}}\)
Esta é a distribuição Beta \(\text{Beta}(1/2, 1/2)\), que é própria porque a integral da densidade de Jeffreys sobre o intervalo \([0, 1]\) é finita.
Temos que a verossimilhança é:
\(L(\theta | x_1, \ldots, x_n) = \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i}\)
Com a priori de Jeffreys, a posteriori resulta em:
\(\pi(\theta | x_1, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i} (1 - \theta)^{n - \sum_{i=1}^n x_i} \cdot \frac{1}{\sqrt{\theta(1 - \theta)}}\)
\(\pi(\theta | x_1, \ldots, x_n) \propto \theta^{\sum_{i=1}^n x_i - 1/2} (1 - \theta)^{n - \sum_{i=1}^n x_i - 1/2}\)
Esta é uma distribuição Beta \(\text{Beta}(\sum_{i=1}^n x_i + 1/2, n - \sum_{i=1}^n x_i + 1/2)\).
O estimador bayesiano sob perda quadrática é a média da distribuição a posteriori Beta:
\(\hat{\theta}_{\text{PQ}} = \mathbb{E}[\theta | x_1, \ldots, x_n] = \frac{\sum_{i=1}^n x_i + 1/2}{n + 1}\)
Para encontrar o intervalo HPD, precisamos calcular os quantis da distribuição Beta posteriori \(\text{Beta}(\sum_{i=1}^n x_i + 1/2, n - \sum_{i=1}^n x_i + 1/2)\).
Como a questão não explicita o total de observações e a soma das observações, podemos visualizar o comportamento do HPD ao passo em que o tamanho amostral aumenta:
def hpd_interval(a, b, level=0.99, sample_size=100000):
sample = beta.rvs(a, b, size=sample_size)
sample_sorted = np.sort(sample)
interval_idx_inc = int(np.floor(level * sample_size))
n_intervals = sample_size - interval_idx_inc
interval_width = sample_sorted[interval_idx_inc:] - sample_sorted[:n_intervals]
min_idx = np.argmin(interval_width)
hpd_min = sample_sorted[min_idx]
hpd_max = sample_sorted[min_idx + interval_idx_inc]
return hpd_min, hpd_max
def update_plot(n):
media_x = 0.7 # Supondo uma média dos xi
sum_x_i = media_x * n
a = sum_x_i + 0.5
b = n - sum_x_i + 0.5
x = np.linspace(0, 1, 1000)
y = beta.pdf(x, a, b)
hpd_99 = hpd_interval(a, b, level=0.99)
hpd_99_formatted = (round(hpd_99[0], 3), round(hpd_99[1], 3))
plt.figure(figsize=(6, 3))
plt.plot(x, y, label=f'Beta({a:.2f}, {b:.2f})')
plt.fill_between(x, y, where=(x >= hpd_99[0]) & (x <= hpd_99[1]), color='gray', alpha=0.5, label='HPD')
plt.title(f'Distribuição Beta e Intervalo HPD de 99%\nIntervalo HPD: {hpd_99_formatted}')
plt.xlabel('θ')
plt.ylabel('Densidade de Probabilidade')
plt.legend()
plt.show()
n_slider = IntSlider(min=1, max=100, step=1, value=10, description='n')
interactive_plot = interactive(update_plot, n=n_slider)
display(interactive_plot)O gráfico supracitado foi feito pressupondo que a média de \(x_i\) fosse 0.7 , pode-se perceber que a média da distribuição posteriori se torna mais estável e menos sensível a pequenas alterações nos dados conforme 𝑛 aumenta. Isso significa que a estimativa pontual de 𝜃 se torna mais confiável.
Os parâmetros \(a\) e \(b\) da distribuição Beta aumentam linearmente com \(n\), tornando a distribuição mais estreita. Para tamanhos amostrais pequenos, a distribuição Beta toma formas largas e até assimétricas, dependendo dos valores observados. Conforme \(n\) aumenta, a distribuição tende a se tornar mais simétrica e mais próxima de uma distribuição normal centrada na média posteriori e ,consequentemente, reduzindo o comprimento do intervalo HPD, convergindo para um valor fixo muito próximo da média verdadeira de \(\theta\).
Considere uma amostra \(y_1, \ldots, y_n\) da distribuição Normal com média \(\mu\) e variância \(\sigma^2 = 1/\tau\) desconhecidas, e suponha que a priori a distribuição de \((\mu, \tau)\) é especificada da seguinte forma: \(\mu | \tau \sim N(\mu_0, 1/\lambda_0 \tau)\) e \(\tau \sim \text{Gama}(\alpha_0, \beta_0)\), onde \(\lambda_0\), \(\alpha_0\) e \(\beta_0\) são positivas.
Ressaltando:
A função de verossimilhança para uma amostra \(y_1, \ldots, y_n\) é:
\(p(y | \mu, \tau) = \left( \frac{\tau}{2\pi} \right)^{n/2} \exp\left( -\frac{\tau}{2} \sum_{i=1}^n (y_i - \mu)^2 \right)\)
Temos a posteriori:
\(p(\mu, \tau | y) \propto p(y | \mu, \tau) p(\mu | \tau) p(\tau)\)
Substituindo as expressões:
\[ \begin{align*} p(\mu, \tau \mid y) &\propto \left( \frac{\tau}{2\pi} \right)^{n/2} \exp\left( -\frac{\tau}{2} \sum_{i=1}^n (y_i - \mu)^2 \right) \\ &\quad \times \left( \frac{\lambda_0 \tau}{2\pi} \right)^{1/2} \exp\left( -\frac{\lambda_0 \tau}{2} (\mu - \mu_0)^2 \right) \\ &\quad \times \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)} \tau^{\alpha_0 - 1} \exp(-\beta_0 \tau) \end{align*} \]
Agrupando os termos dependentes de \(\mu\) e \(\tau\):
\(p(\mu, \tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp\left( -\frac{\tau}{2} \left[ \sum_{i=1}^n (y_i - \mu)^2 + \lambda_0 (\mu - \mu_0)^2 \right] - \beta_0 \tau \right)\)
Para integrar \(\mu\), completamos o quadrado no expoente:
\(\sum_{i=1}^n (y_i - \mu)^2 + \lambda_0 (\mu - \mu_0)^2\)
Expandindo e agrupando os termos:
\(\sum_{i=1}^n y_i^2 - 2\mu \sum_{i=1}^n y_i + n\mu^2 + \lambda_0 \mu_0^2 - 2\lambda_0 \mu \mu_0 + \lambda_0 \mu^2\)
Agrupando os termos quadráticos, lineares e constantes em \(\mu\):
\(\left( n + \lambda_0 \right) \mu^2 - 2 \left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right) \mu + \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2\)
Completando o quadrado:
\(\left( n + \lambda_0 \right) \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 + \left( \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 \right) - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0}\)
Assim, o expoente se torna:
\(-\frac{\tau}{2} \left[ (n + \lambda_0) \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 + \sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0} \right]\)
A integral em \(\mu\) da expressão acima é uma integral gaussiana, logo:
\(\int_{-\infty}^{\infty} \exp \left( -\frac{\tau (n + \lambda_0)}{2} \left( \mu - \frac{\sum_{i=1}^n y_i + \lambda_0 \mu_0}{n + \lambda_0} \right)^2 \right) d\mu = \sqrt{\frac{2\pi}{\tau (n + \lambda_0)}}\)
Portanto, ao desprezarmos essa constante, a distribuição marginal de \(\tau\) é:
\(p(\tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp \left( -\tau \left( \frac{\sum_{i=1}^n y_i^2 + \lambda_0 \mu_0^2 - \frac{\left( \sum_{i=1}^n y_i + \lambda_0 \mu_0 \right)^2}{n + \lambda_0}}{2} + \beta_0 \right) \right)\)
Simplificando a constante no expoente:
\(p(\tau | y) \propto \tau^{(n+1)/2 + \alpha_0 - 1} \exp \left( -\tau \left( \frac{\sum_{i=1}^n (y_i - \bar{y})^2}{2} + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{2(\lambda_0 + n)} + \beta_0 \right) \right)\)
Reconhecendo a forma da distribuição Gama, temos:
\(\tau | y \sim \text{Gama} \left( \alpha_0 + \frac{n}{2}, \beta_0 + \frac{1}{2} \left[ \sum_{i=1}^n (y_i - \bar{y})^2 + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{\lambda_0 + n} \right] \right)\)
Dada a \(\tau\), a distribuição condicional de \(\mu\) é:
\(\mu | \tau, y \sim N \left( \frac{\lambda_0 \mu_0 + n \bar{y}}{\lambda_0 + n}, \frac{1}{\tau (\lambda_0 + n)} \right)\)
\(\tau | y \sim \text{Gama} \left( \alpha_0 + \frac{n}{2}, \beta_0 + \frac{1}{2} \left[ \sum_{i=1}^n (y_i - \bar{y})^2 + \frac{\lambda_0 n (\bar{y} - \mu_0)^2}{\lambda_0 + n} \right] \right)\)
Dada a \(\tau\), a distribuição condicional de \(\mu\) é:
\(\mu | \tau, y \sim N \left( \frac{\lambda_0 \mu_0 + n \bar{y}}{\lambda_0 + n}, \frac{1}{\tau (\lambda_0 + n)} \right)\)
Essa priori é escolhida para ser não informativa, significando que temos pouca ou nenhuma informação a priori sobre \(\mu\) e \(\tau\).
Distribuição a posteriori: Com a priori não informativa, a posteriori ainda é própria, pois as integrais de normalização convergem.
Comparação com Resultados Clássicos:
A distribuição marginal de \(\mu\) se aproxima da distribuição normal clássica para a média de uma amostra, especialmente à medida que \(n\) aumenta.
A distribuição marginal de \(\tau\) se aproxima da distribuição qui-quadrado clássica (quando reescalada), que é usada para inferência sobre a variância na estatística frequentista.
Seja \(n = (n_1, \ldots, n_k)\) um vetor aleatório com distribuição multinomial e densidade \(p(n | \theta) \propto \prod_{i=1}^k \theta_i^{n_i}\), onde \(\theta = (\theta_1, \ldots, \theta_k)\), \(\theta_i > 0\) e \(\sum_{i=1}^k \theta_i = 1\). Considere a priori para \(\theta\) uma distribuição de Dirichlet com parâmetro \(\alpha = (\alpha_1, \ldots, \alpha_k)\), isto é \(p(\theta) \propto \prod_{i=1}^k \theta_i^{\alpha_i - 1}\).
No caso particular \(\alpha = (0, 0, \ldots, 0)\) a distribuição a priori de \(\theta\) é imprópria. Mostre que a distribuição a posteriori é própria se e somente se \(n_i > 0\) para \(i = 1, 2, \ldots, k\).
Quando \(\alpha = (0, 0, \ldots, 0)\), a distribuição a priori de \(\theta\) é:
\(p(\theta) \propto \prod_{i=1}^k \theta_i^{-1}\)
Esta distribuição é imprópria porque a integral sobre o simplex \(\{\theta \in \mathbb{R}^k : \theta_i > 0, \sum_{i=1}^k \theta_i = 1\}\) não converge.
A função de verossimilhança da distribuição multinomial é:
\(p(n | \theta) \propto \prod_{i=1}^k \theta_i^{n_i}\)
Temos assim a posteriori:
\(p(\theta | n) \propto p(n | \theta) p(\theta) \propto \left( \prod_{i=1}^k \theta_i^{n_i} \right) \left( \prod_{i=1}^k \theta_i^{-1} \right) = \prod_{i=1}^k \theta_i^{n_i - 1}\)
Esta é a forma da distribuição Dirichlet com parâmetros \((n_1, n_2, \ldots, n_k)\):
\(p(\theta | n) \propto \prod_{i=1}^k \theta_i^{n_i - 1}\)
Para que a distribuição a posteriori seja própria, a integral da densidade da posteriori sobre o simplex deve ser finita. A distribuição Dirichlet \(\text{Dirichlet}(n_1, n_2, \ldots, n_k)\) é própria se e somente se todos os \(n_i > 0\). Se qualquer \(n_i = 0\), o termo \(\theta_i^{-1}\) aparece na densidade, e a integral não é finita.
Portanto, a distribuição a posteriori é própria se e somente se \(n_i > 0\) para todos os \(i = 1, 2, \ldots, k\).
Na véspera do primeiro turno para a eleição de governador do DF de 2010, a Datafolha divulgou uma pesquisa indicando que, de 891 eleitores entrevistados que já tinham decidido em quem votar, Agnelo Queiroz tinha a preferência de 467, Weslian Roriz a de 315 e outros candidatos a de 109 eleitores.
Formule um modelo para analisar esses dados. O interesse centra fundamentalmente em três perguntas:
Usaremos um modelo baseado na distribuição Multinomial para representar os votos dos eleitores. A distribuição Multinomial é adequada aqui porque modela a probabilidade de um certo número de sucessos em várias categorias.
Seja \(\theta = (\theta_1, \theta_2, \theta_3)\) as proporções de votos para Agnelo, Weslian e outros candidatos, respectivamente, onde:
\(\theta_1 + \theta_2 + \theta_3 = 1\)
Os dados \(n = (n_1, n_2, n_3)\) seguem uma distribuição multinomial:
\((n_1, n_2, n_3) \sim \text{Multinomial}(n, \theta)\)
Para a análise bayesiana, assumimos uma distribuição a priori para \(\theta\). Podemos usar uma distribuição de Dirichlet com parâmetros \(\alpha = (\alpha_1, \alpha_2, \alpha_3)\):
\(\theta \sim \text{Dirichlet}(\alpha)\)
Dado que os dados seguem uma distribuição multinomial e a priori é uma Dirichlet, a distribuição a posteriori de \(\theta\) também é uma Dirichlet com parâmetros \(\alpha' = (\alpha_1 + n_1, \alpha_2 + n_2, \alpha_3 + n_3)\):
\(\theta | n \sim \text{Dirichlet}(\alpha_1 + n_1, \alpha_2 + n_2, \alpha_3 + n_3)\)
Para uma priori não informativa, podemos assumir \(\alpha_i = 1\) para todos os \(i\).
Para que a eleição seja definida no primeiro turno, um candidato deve obter mais de 50% dos votos válidos. Vamos calcular a probabilidade posteriori de \(\theta_1 > 0.5\):
\(P(\theta_1 > 0.5 | n)\)
Semelhante à (a), precisamos calcular:
\(P(\theta_1 > 0.5 | n)\)
Calculamos a diferença na proporção dos votos válidos entre Agnelo e Weslian:
\(P(|\theta_1 - \theta_2| > d | n)\)
onde \(d\) é a diferença que desejamos calcular.
Podemos implementar isso usando simulações a partir da distribuição posteriori Dirichlet para estimar essas probabilidades.
n_total = 891
n_agnelo = 467
n_weslian = 315
n_outros = 109
alpha = np.array([1, 1, 1])
alpha_post = alpha + np.array([n_agnelo, n_weslian, n_outros])
n_samples = 100000
samples = dirichlet(alpha_post).rvs(n_samples)
prob_agnelo_eleito = np.mean(samples[:, 0] > 0.5)
dif_agnelo_weslian = samples[:, 0] - samples[:, 1]
diff_percent = np.mean(dif_agnelo_weslian)
print(f"Probabilidade de Agnelo ser eleito no primeiro turno: {prob_agnelo_eleito:.2f}")
print(f"Diferença média na porcentagem entre Agnelo e Weslian: {diff_percent:.2f}")
plt.figure(figsize=(5,5))
x = np.linspace(0, 1, 1000)
plt.plot(x, beta.pdf(x, alpha_post[0], alpha_post[1] + alpha_post[2]), label="Agnelo")
plt.plot(x, beta.pdf(x, alpha_post[1], alpha_post[0] + alpha_post[2]), label="Weslian")
plt.axvline(0.5, color='red', linestyle='--', label="50% Threshold")
plt.xlabel("Proporção")
plt.ylabel("Densidade de Probabilidade")
plt.title("Distribuições Posteriori das Proporções de Votos")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper right')
plt.grid(True)
plt.tight_layout()
plt.show()Probabilidade de Agnelo ser eleito no primeiro turno: 0.92
Diferença média na porcentagem entre Agnelo e Weslian: 0.17
Com base nos dados da pesquisa e na análise bayesiana, é altamente provável que Agnelo Queiroz seja eleito no primeiro turno, com uma vantagem média de 17% sobre Weslian Roriz. A probabilidade de Agnelo obter mais de 50% dos votos válidos é muito alta (91.97%), o que sugere que a eleição pode ser decidida no primeiro turno.
Considere \(n = 12\) ensaios de Bernoulli com probabilidade de sucesso \(p\) e \(y = 9\) sucessos. Suponha que a distribuição a priori de \(p\) é especificada de forma que \(\eta = \log \left( \frac{p}{1 - p} \right)\) segue uma distribuição Normal com média \(\mu = 0\) e variância \(\sigma^2 = 1\).
De imediato temos:
\(\eta \sim N(0, \sigma^2)\)
O parâmetro \(p\) pode ser expresso em termos de \(\eta\) como:
\(p = \frac{e^\eta}{1 + e^\eta}\)
A verossimilhança para a distribuição Bernoulli é:
\(p(y | p) = \binom{n}{y} p^y (1 - p)^{n - y}\)
Para \(n = 12\) e \(y = 9\):
\(p(9 | p) = \binom{12}{9} p^9 (1 - p)^3\)
Para calcular \(\mathbb{E}(p | y = 9)\) e \(\text{DP}(p | y = 9)\), podemos usar a aproximação de Laplace, que envolve encontrar a moda da posteriori e aproximar a distribuição posteriori por uma Normal em torno da moda.
Podemos implementar isso em Python para calcular as aproximações solicitadas.
n = 12
y = 9
mu = 0
variances = [1, 4, 9]
def posterior_approximation(y, n, mu, sigma2):
# Log-verossimilhança em termos de eta
def log_likelihood_eta(eta):
p = expit(eta)
return y * np.log(p) + (n - y) * np.log(1 - p)
def prior_eta(eta):
return norm.logpdf(eta, loc=mu, scale=np.sqrt(sigma2))
def log_posterior_eta(eta):
return log_likelihood_eta(eta) + prior_eta(eta)
def dlog_posterior_eta(eta):
p = expit(eta)
return y - n * p - eta / sigma2
def d2log_posterior_eta(eta):
p = expit(eta)
return -n * p * (1 - p) - 1 / sigma2
# Modo da posteriori (máximo a posteriori)
from scipy.optimize import minimize
result = minimize(lambda eta: -log_posterior_eta(eta), x0=0, method='Newton-CG', jac=lambda eta: -dlog_posterior_eta(eta), hess=lambda eta: -d2log_posterior_eta(eta))
eta_mode = result.x[0]
# Variância em torno da moda
eta_var = -1 / d2log_posterior_eta(eta_mode)
eta_samples = np.random.normal(loc=eta_mode, scale=np.sqrt(eta_var), size=10000)
p_samples = expit(eta_samples)
mean_p = np.mean(p_samples)
std_p = np.std(p_samples)
prob_p_greater_half = np.mean(p_samples > 0.5)
return p_samples, mean_p, std_p, prob_p_greater_half
fig, axes = plt.subplots(len(variances), 1, figsize=(5,10))
for ax, sigma2 in zip(axes, variances):
p_samples, mean_p, std_p, prob_p_greater_half = posterior_approximation(y, n, mu, sigma2)
ax.hist(p_samples, bins=50, density=True, alpha=0.75)
ax.axvline(mean_p, color='r', linestyle='--', linewidth=2, label=f'Média = {mean_p:.3f}')
ax.axvline(0.5, color='g', linestyle='-.', linewidth=2, label='$p = 0.5$')
ax.set_title(f'Posterior de $p$ para $\sigma^2 = {sigma2}$', fontsize=14)
ax.set_xlabel('$p$', fontsize=12)
ax.set_ylabel('Densidade', fontsize=12)
ax.legend(loc='upper left', fontsize=12)
ax.grid(True)
print(f"Para σ² = {sigma2}:")
print(f" E(p | y = 9) = {mean_p:.4f}")
print(f" DP(p | y = 9) = {std_p:.4f}")
print(f" Pr(p > 0.5 | y = 9) = {prob_p_greater_half:.4f}")
print()
plt.tight_layout()
plt.show()Para σ² = 1:
E(p | y = 9) = 0.6736
DP(p | y = 9) = 0.1099
Pr(p > 0.5 | y = 9) = 0.9275
Para σ² = 4:
E(p | y = 9) = 0.7144
DP(p | y = 9) = 0.1168
Pr(p > 0.5 | y = 9) = 0.9502
Para σ² = 9:
E(p | y = 9) = 0.7240
DP(p | y = 9) = 0.1208
Pr(p > 0.5 | y = 9) = 0.9469
Percebemos que a medida que a variância da priori aumenta, a incerteza na posteriori também aumenta, refletindo uma maior incerteza a priori sobre \(p\), sendo propagada para a distribuição a posteriori.